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It has been well-known that under the assumption of a constant uniform mean 
flow, the acoustic wave propagation equation can be formulated as a boundary inte- 
gral equation, in both the time domain and the frequency domain. Compared with 
solving partial differential equations, numerical methods based on the boundary 
integral equation have the advantage of a reduced spatial dimension and, hence, re- 
quiring only a surface mesh. However, the constant uniform mean flow assumption, 
while convenient for formulating the integral equation, does not satisfy the solid 
wall boundary condition wherever the body surface is not aligned with the uniform 
mean flow. In this paper, we argue that the proper boundary condition for the 
acoustic wave should not have its normal velocity be zero everywhere on the solid 
surfaces, as has been applied in the literature. A careful study of the acoustic energy 
conservation equation is presented that shows such a boundary condition in fact 
leads to erroneous source or sink points on solid surfaces not aligned with the mean 
flow. A new solid wall boundary condition is proposed that conserves the acoustic 
energy and a new time domain boundary integral equation is derived. In addition 
to conserving the acoustic energy, another significant advantage of the new equation 
is that it is considerably simpler than previous formulations. In particular, tangen- 
tial derivatives of the solution on the solid surfaces are no longer needed in the 
new formulation, which greatly simplifies numerical implementation. Furthermore, 
stabilization of the new integral equation by Burton-Miller type reformulation is 
presented. The stability of the new formulation is studied theoretically as well as 
numerically by an eigenvalue analysis. Numerical solutions are also presented that 
demonstrate the stability of the new formulation. 


I. Introduction 


The numerical solution of sound scattering by an acoustically large body remains a significant 
challenge due to the high demand on computational resources that are required to resolve the 
acoustic waves of short wavelengths. It is well-known that under the assumption of a constant 
mean flow, the acoustic wave propagation is governed by the convective wave equation that, in turn, 
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can be converted into a boundary integral equation. The boundary integral equation approach has 
the advantage of reducing the spatial dimensions of the problem by one, making it an attractive 
computational method for calculating sound scattering and shielding at mid to high frequencies. 
In this paper, we consider the problem of acoustic scattering by rigid bodies in the presence of a 
uniform flow using the boundary integral equation approach. The present approach is based on the 
time domain boundary integral equation. The time domain approach has some distinct advantages 
over a frequency domain approach. Most notably, scattering solutions at all frequencies are obtained 
within one single computation. In addition, broadband noise sources and time dependent transient 
signals can be simulated and studied. The time domain approach also couples naturally with 
nonlinear computations where many frequencies are generated. 


Previously, scattering of sound waves by rigid bodies with flow has been studied, in both the 
frequency domain and the time domain. In [30], acoustic radiation in moving flow was formulated 
as a boundary integral equation in the frequency domain. The nonuniqueness of the exterior problem 
was dealt with by applying the Burton-Miller reformulation procedure.! The time domain boundary 
integral equation approach for scattering by moving surfaces was first formulated and studied in 
[28]. More recent studies of the time domain approach in the presence of a mean flow can be found 
in [13, 16, 19]. 

A major difference between the current approach and those taken previously is in the treatment of 
the boundary condition at solid surfaces in the presence of flow. While the linear acoustic problem 
as a perturbation over the mean flow can be considered separately from the mean flow, an implicit 
condition is that the mean flow itself satisfies the solid wall boundary condition. The assumption of 
a constant mean flow is an approximation to the actual mean flow and this assumption is made to 
make the formulation of a boundary integral equation possible. While this facilitates the conversion 
of the partial differential equation to the boundary integral equation, the simplified mean flow itself 
obviously can not satisfy the physical boundary condition at solid boundaries wherever the surface 
is not aligned with the assumed constant mean flow. As pointed out in [28], the boundary integral 
equation derived based on such an assumption would be formally valid when M,, < 1 where M,, is 
the Mach number of mean flow normal to the body surface. In this paper, we take a closer look at 
the boundary condition to be used for scattering of acoustic waves at solid surfaces where M,, is not 
zero. In all the previous studies, a boundary condition of normal acoustic velocity being zero has 
been applied everywhere including the surfaces where M,, 4 0. However, an analysis of the acoustic 
energy equation shows that the usual boundary condition would lead to nonzero energy flux at 
surfaces where M, # 0. A new formulation is derived based on this acoustic energy consideration, 
and an alternative boundary condition is proposed by the requirement that energy flux be zero at 
solid surfaces. As we will see, due to the structure of the integral equation, the new formulation 
also becomes much simpler than those found in the literature for scattering with flow, which is of 
great benefit for computation. 


In addition to the modification of boundary condition at solid surfaces, a Burton-Miller type re- 
formulation of the integral equation consistent with the new boundary condition is also presented. 
It is well-known that the direct solution of the boundary integral equation for exterior scattering 
problems is prone to numerical instabilities. 2% !%!%18:30 In the time domain, the instability is 
more easily excited because all frequencies within the numerical resolution are present. There are 
generally two approaches for dealing with this instability. One is the Burton-Miller reformulation 
which has been widely used for frequency domain exterior scattering problems. Recently, it has 
been shown that Burton-Miller reformulation is effective for time domain as well.2%!8 Another 
method for the removal of the instability is the CHIEF method.!"?9 In the present study, we apply 
the Burton-Miller technique to our new formulation for the elimination of instabilities. 
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The rest of the paper is organized as follows. In Section 2, an integral relation for acoustic wave 
propagation is derived for a constant mean flow in a general direction. Then, the time domain 
boundary integral equation for scattering by rigid bodies is derived in Section 3. In Section 4, a 
Burton-Miller type reformulation of time domain boundary integral equation is presented and a 
discussion on the stability of the new formulation is given in Section 5. Numerical methods for 
the time domain boundary integral equation are discussed in Section 6. Stability of the current 
formulation is demonstrated in Section 7 by analyzing the eigenvalues of the discretized system. 
Section 8 contains the conclusions. 


II. Integral representation of acoustic waves in the presence of a uniform 
mean flow 


The current problem is considered in the context of solving the wave equation in a moving medium 
exterior of certain specified surface S, such as the scattering of sound field by an object as shown 
in Figure 1. Acoustic waves are assumed to be disturbances of small amplitudes. Linear acoustic 
problems are frequently formulated using a velocity potential function ¢(r,t) where the acoustic 
velocity u and pressure p are related to ¢ as follows: 


u=Vé, p=-m(Z+u-ve) (1) 


where po is the mean density. With a constant mean flow U, the acoustic disturbances are governed 
by the convective wave equation.”° In the present study, we consider the solution of the following 
equation for the velocity potential: 


2 
(= Fe ce v) 6 — 2V9 = (r,t) (2) 
with homogeneous initial conditions 
o(r, 0) = 2 (n.0) =0, t=0 (3) 


In the above, c is the speed of sound, U is the constant mean velocity, and q(r,t) represents the 
known acoustic sources. Furthermore, in addition to the radiation condition at the far field, (2) 
and (3) are to be supplemented with boundary conditions on the scattering surface S. The suitable 
boundary conditions to be applied on solid surfaces will be discussed in Section 3. 


Sources Ss 


ae 


Figure 1. A schematic showing the scattering body and mean flow. Scattering surface is denoted by 
S and the solution domain exterior of S is denoted by V. The surface normal vector n is taken to be 
inward toward the interior of the body. 
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It is well-known that the convective wave equation (2) and the initial condition (3) as well as the 
boundary conditions can be reformulated into an integral equation. In the literature, the integral 
representation of sound waves in a moving flow is often derived by making use of generalized 
functions in a setting of moving bodies in an otherwise undisturbed medium.*%*?0:2%:26 Here, 
we present a derivation using a free-space Green’s function G(r, t; r’,t') that, for convenience of 


discussion, is defined as follows: 


2 
(5 +u-v) G-—V’G = 6(r—r')5(t —t’) (4) 
with initial conditions 
G(r .tsr!,t) =F (w,ts9!,t) =0, GST. (5) 


Note that the time domain Green’s function G(r, t;r’, t/) defined above is nonzero for t € (—oo, t’]. 
The solution to (4) and (5) is well-known (see, e.g., [4, 10, 25]) and, for a mean flow of a general 
direction, can be written as 


. Go R 
G(rstin't) = 2% (¢ +8: (1) - 5) (6) 
where 
G : dR = \/[M \P +02 Me 7 
0= Fey amd Reve!) = VIM rr? + 02 |n — ot (7) 
in which 
U U U M 
2, = VI-M, B= =a a= URL, M=IM| 8) 


By an operation of Gx (2) —¢x(4) and by integrating over the volume V exterior of the scattering 
surface S for space and an interval [0~,¢’*] for time ¢, it is straight-forward to show that we will 
get 


(a(+u-ve) -6( FZ av-va))u 


_ ey. |eve 2 ov \ drdt = i. a |Gar, t) — o(r, t)d(r — r')5(t — | drdt 


Integration of the first term in the above will be zero by initial conditions thus defined for ¢ and 


G. Then, upon using the divergence theorem and the condition at infinity, we get an expression for 
pressure @ at an arbitrary point r’ in V and time ?’ as follows: 
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ut tt 


o(r ‘=f [Gate narar+e [ [ese -05 2G ar sat 
ik 6 (Bru. v6) - (Gi +u ve) 


where r, denotes points on surface S$, and 


M,,dr ,dt (9) 


M,=n:-M=n-U/c 


is the normal component of the mean velocity Mach number on surface point r,. Here, the unit 
normal vector n is assumed to be outward from the solution domain. For the exterior scattering 
problem considered in the present study, the normal vector is then the one that is inward to the 
body as noted in Figure 1. 


For convenience of discussion, we define a modified normal derivative (denoted by an overbar) as 


a a 
an = 9q 7 Ma(M-Y) (10) 


Then, equation (9) can be written as 


t+ es 
wlth = [J Gatesrararne [ [G20 rdrsare [" [ [652 — 02] anarat 
(11) 
Furthermore, if we introduce a combined normal derivative (denoted by a tilde) as 
0 0 6M, [0 0 6M, 0 
= : sce ee, Moe 12 
On Bn (3 none v) = One Ot a) 


we get another expression: 


ot) = ff Garnara+e f 


Equation (9), (11) or (13) is the Kirchhoff integral representation of the acoustic field in the presence 
of a uniform mean flow. The integral relation can be further expressed as the integration of retarded 
values by utilizing G as given in (6). In particular, note that we have 


tit tt 


[eg - 95 (Oar sil (13) 


[5 (e248 r) a) + a (¢ t+B-(r'-7)-=5)| (14) 


where Go and R are those defined in (7). Then equation (13) can be written as 


OG 1 AG 
On Arc? On 


ry 1 1 0 OG ; R 
Blt) = ate ff patrstndans ff [GoSE (rate) — 3 (oleate) + eS (ra te) dr 


(15) 
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where V, denotes the region of acoustic sources and the retarded time for t’ is defined as 


R 


teat +B-(r—r)- 5 eo 
The modified normal derivative for Go is found to be the following: 
OGo 1 OR on-(r—r’) 
— =a — = 1 
on BeOn - R8 oe 


Equation (15) relates the solution at point r’ and time ?t’ to the direct contribution from source 
function q and a surface contribution involving the retarded values of ¢ and their normal derivatives. 
As shown in [13], this form is equivalent to previous such formulations appearing in the literature, 
e.g., in [23,26], where the relationship had been derived under the assumption of a mean flow that 
is aligned with the z-axis. 


When both ¢(rz,t) and 20 (pt) on surface S are known, ¢(r’,t’) at any field point r’ can be 
computed by using (15). However, (r,t) and 28 (r, t) are not independent. They have to satisfy 


the boundary integral equation formed when r’ is taken to be a boundary point r‘, as we will discuss 
next. 


III. Time domain boundary integral equation for scattering with solid surfaces 


A Boundary Integral Equation (BIE) is formed by taking the limit r’ > r/, in the integral relation 
(15), where r‘, is a point on the boundary. The integral in (15) involving Go is weakly-singular 
and, by using equation (55) given in the Appendix (assuming r’, is a smooth boundary collocation 
point), it can be shown that 


tim, [Sree olrosteddrs= f Era )6(rsta)drs—26(r4,0) (18) 
S 


rr, Sg on 
Applying this limit to (15), we get the following Time Domain Boundary Integral Equation (TD- 
BIE): 


! ! 0 ! OG / R 0 / / / 
2notrist)— f (GoFE Crate) — F2 [orate + a Blrat)] ) dr = Ort) 9) 


where Q(r%, t’) denotes the contribution from the external sources to the surface point r/,: 


1 1 
1 ahh ! 
Qe =f palrtadar (20) 
For sound scattering problems, ¢(r%, t’) on the scattering surface S is to be determined by (19) when 
the boundary condition for ¢ on S is given. A customary boundary condition on rigid surfaces is 
that the normal component of the acoustic velocity is zero, i.e., w-n = 0, which, considering (1), 
leads to 


n-Vo= A(r,,t) =0, rseS (21) 


Indeed, in all the previous literature on wave scattering with a uniform mean flow (e.g., [5, 10, 
13, 15, 19, 28, 30]), in both the frequency domain and the time domain, boundary conditions of 
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type (21) have been assumed at solid wall boundaries. To implement such a boundary condition, 
the combined normal derivative appearing in (19) would then be separated into the normal and 
tangential components as 


Od 9 Ob 10¢ 
=(1-—M, M, +} Mr-V 22 
On ( n) On "\c at rt Ve ) 
where M7 is the tangential component of the mean flow Mach number M. 

In the present paper, however, we propose an alternative boundary condition to be used at solid 
surfaces when solving TDBIE (19) in the presence of a uniform flow. The new boundary condition 
is based on a consideration of the acoustic energy. 


It can be shown that the convective wave equation (2) without the source term has an associated 
energy equation: 


JE 
ae . —*, 2 
atv Jy =0 (23) 
where 
low. 1 [De U-VdD¢ _ 06 1 Dé D_ oO 
ple + 92 | Di 2 pr ae ep pe Ee 2s) 


Equation (23) can be validated directly by using the expressions defined in (24). When substituted 
by the acoustic velocity and pressure defined in (1), po is the usual acoustic energy density in a 
uniform flow.22 24:27 


By (24), it is immediately clear that the energy flux at a surface of normal n is the following: 


2 _ 0¢ (06 M,De\ _ d¢0¢ 
oi Ma (3 C a) = Ot Oi oo) 


Clearly, on a surface where the normal component of the mean velocity M), is not zero, i.e., where 
the surface is not aligned with the mean flow, application of boundary condition (21) will result in 
nonzero energy flux, i.e., J, 4 0 and, consequently, cause the surface to be acting like an acoustic 
energy source or sink according to (25). This will apparently lead to nonconservation of the total 
acoustic energy. 


Alternatively, the boundary condition on the solid surface may be defined by the requirement that 
no energy flows into or out of the surface. By (25) and to ensure energy flux J; = 0 on solid surfaces, 
we propose that the boundary condition be modified such that the combined normal derivative of 
@, defined in (12), is zero: 


0¢ M,D¢ 

= t) = — - —— =0 S 26 
(rat) = So - ATE =0, ree (26) 
The total acoustic energy will be conserved under this new condition. 

Now by applying boundary condition (26) to (19), a new formulation of the TDBIE for ¢(r%, t’) 


with solid surfaces is found as follows: 


OG RO 
Qnb(r’,,t') + [ a (ortho + arate) dr, = Q(r",, t’) (27) 
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Equation (27) is one of the main results of the present paper. It is a new formulation for the time 
domain boundary integral equation for acoustic scattering by rigid surfaces in a constant mean flow. 
It is different from those in the literature in several aspects. First, the boundary condition used for 
(27) is one that is based on the acoustic energy flux consideration instead of the acoustic normal 
velocity. The two approaches differ on the part of the boundary where the mean flow itself does 
not satisfy the slip boundary condition. Second, the new equation is much simpler than those of 
the previous formulations in which tangential derivatives of the solution on the scattering surface 
are required to be kept as part of the integral equation. Of course, boundary condition (26) reduces 
to the usual one (21) wherever the mean flow does satisfy the solid wall boundary condition, i.e., 
My, = 0. 


IV. Burton-Miller type reformulation in time domain with a mean flow 


Direct solution of boundary integral equations for exterior scattering problems, however, is known 
to suffer numerical instabilities. The instability is generally attributed to the existence of resonance 
frequencies for the interior domain.!*:% 3° In time domain solutions, the instability is more easily 
triggered as a continuous spectrum of frequencies within the numerical resolution are present in the 
computation. This instability is one of the major difficulties that have hindered the use of time 
domain integral equations. Recently, Burton-Miller type reformulation that has been widely used 
for exterior scattering problems in the frequency domain has shown to be effective in eliminating the 
instability in the time domain as well.!:®?° In [2], a theoretical justification has been provided for 
the extension of the Burton-Miller formulation to the time domain for the wave equation without 
flow. In this section, we propose the Burton-Miller reformulation for the TDBIE (27). An analysis 
on its stability similar to that in [2] is given in the next section. 


For convenience of discussion, we define the following time domain double layer potential: 
t+ 


dG Ree _ f dGo ; f R 0¢ ; 
[ Gerster trotrs tarde = [Fer (orate) + Crate) drs 


s 
(28) 


The Burton-Miller type reformulation is carried out by applying a linear combination of the time 


Diet.) = f 


and certain normal derivatives to the time domain integral equation. In earlier studies of the 
Burton-Miller formulation for scattering with a flow, the modified normal derivative (10) had been 
used.!*:3° Here, we propose that the normal derivative to be used for the Burton-Miller formulation 
be the combined normal derivative defined in (17). Specifically, the Burton-Miller reformulation is 
obtained by applying the following derivative operator to the boundary integral equation at surface 
points 7: 


fhe 2 
a= + bc A (29) 


where @ and 6 are constants and c is the speed of sound, namely, 


2 (4no(o',#) + Dio\r',1)) 


On! 


0 / / / / re 
a=, (2r6(r4,t'}+DIdl(,1')) +be 


(30) 
Applying again the solid surface boundary condition (26), equation (30) is expanded to be the 
following: 
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SPB Sen DE pf OR.. tihe c REO 
a lan Z (wt) +f Seren!) (Berasth) + Ay Solraste)) dr 


pe O OG , ! R Od F _ _O paar - Oo 5 hi 
clam f Bint (Hast? ) (Olrata) + Foret dr. = ae (rst) + bea (rst) 


(31) 


Note that an integral with a kernel 2Go (pr) is hyper-singular when r, coincides with r,’. In 
particular, we have 


0?Go ry O 9n-(rs—1r,) _ az j a[n- (rs —1%,)] [n’ - (vr, —r)| 
On’ On (rs, 7s) = On! | a Re | = Fe [n “nT MyM, + 3a | (32) 


/ 


Thus, 2Go (pr!) is of order O(1/|r; — 1,|3) as rs > r%,. 
We consider the following regularization process for the hyper-singular integral in (31) that adds 


and subtracts a term involving the value at the collocation point ¢(r%, t’): 


O OG 7 / R Og 
dnl | ¥ aq (Ths) (orth a5 re: Dt Fe (rath) ir, 
0 OG j j RO j 0 OG: j 
= | Petra) (erate) ole.) + 5 rath) dra] +6000 | f Ser dry 


(33) 
The first integral is now integrable by Cauchy Principal Value (Appendix B) and the second integral 
is zero according to (55) given in the Appendix A. Upon carrying out the derivatives inside the first 
integral shown above, we get the following Burton-Miller reformulation of the time domain boundary 
integral equation (BM-TDBIE): 


.O¢(ri,t’) . f OGo (9b, R Oo b -30G) 0G) Od, 
eg ee [ an \ Ot (rate) + 3 Ses pez (ste) ) drs — ca [a on on op etnias 


R ad 


G R OQ 
+be , On/On (orate — o(r4,t') + arent) dr, = 


— (1, ,t) + be (r%, #): - (84) 


5 (p! 
“atl an! 


The proper values for the coefficients @ and 6 will be given in the next section where stability of 
(34) will be discussed. 


V. Stability of the time domain Burton-Miller formulation in the presence of 
a mean flow 


Following closely the work in [2] for the case without flow, we demonstrate in this section that 


the Burton-Miller type reformulation presented in the previous section eliminates the nontrivial 
solutions of the homogeneous integral equation in the case with a flow as well. 
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Suppose that there is a nontrivial solution ¢9(r;,t) to the homogeneous formulation (34). We will 
show in what follows that such a solution is not possible. Consider the double layer potential (28) 
extended to domains both exterior and interior of surface S: 


_ wt, rr’ €V, exterior of S 
D Poghy , ! R Oo ! d —- , , 
[go] (1, t) aac cs on er”) go(ts,tr) + a5, (Ts: tr) r= WO; r =f, on §S 
w , 9’ €V_-, interior of S 


We note that wt and w7 satisfy the homogeneous convective wave equation in the exterior and 
interior domains of S, respectively. It can also be shown that 


lim wt = wo — 2r¢o(r4, t’) (35) 
r'or’, 
lim w7 = wo + 27d0(7r%,, t’) (36) 
rr! 
Owt Ow 
li = ii 
ee Bar ot, OF my 


Equations (35) and (36) can be found by using the limits given in (55) in the Appendix, and equation 
(37) follows after an application of the regularization process (33) to both sides of the equation. 


Now since ¢0(15,t) satisfies the homogeneous Burton-Miller formulation (30), we have, at r’ = r%, 


~ O 
270 + wo) + bc Dal (4rrdo + wt) =0 


s 


<a 
“OH 


Tr 


By the jump conditions (35)-(37) as well as the boundary condition (26), the above yields 


Ow + Ow” 


On the other hand, since w~ satisfies the convective wave equation and by the energy equation (23) 
of the convective wave equation, we have 


0 1 ais pe U -Vw~ Dw _ Ow = dpeX 
af |gven © 22 C2 B|er= fv [% (vu ~ @ Dt v)| ar 


which, with an application of the divergence theorem, becomes 


1 1 
=|Vw |? 4 
f [ive 22 


where V~ represents the volume interior of S$. The minus sign on the right hand side has been 
added due to the fact that the normal derivative used in (39) is still the one that is inward of the 
body surface. Note that, for subsonic flows where |U| < c, the left hand side of (39) is nonnegative: 


Dw 
Dt 


Dw- c 


Dt 


dr dt (39) 


U -Vw” Dw r=-[- Ow” Ow 
ce Dt fo dg OF ON 
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1 
2c? 


Dw 
Dt 


2 = = = 
1 _, | Dw U-Vw Dw 
> 
) + vw | Di | 0 


i 9 
=|V 
3 el Cc Dt — 


U-Vw Dw _1 Iw | 1|Dw7 
Cc Dt 2 Dt 


On the other hand, using (38), the right hand side of (39) will be nonpositive: 


tf Ow? Owe 1 ft f @ |dw-|? 
-{ e air = 5 f [el dr, <0 
0 S Ot on Cc 0 S bc Ot 
provided 
220 (40) 
b 
The above implies that w~ has to be a trivial solution, i.e., w~ = 0 under the condition (40). A 
simple choice for @ and bis a= —b=1. 


VI. Time Domain Boundary Element Method 


In this section, we describe a numerical solution of (34) by the Time Domain Boundary Element 
Method (TDBEM) and demonstrate numerical stability of the new formulation. 

Let surface S be discretized by surface elements E;, 7 = 1,2,..., Ne, where Ne is the total number 
of elements, and the time be discretized by t, = nAt, where At is the time step. The time domain 
numerical solution on the surface can be expanded as 


Ne Ne 
O(rs,t) = >) >) uP e;(rs)dn(t) (41) 
n=0 j=l 
where y;(r;) is the surface basis function for element E; and y(t) is the temporal basis function 
for time node t,. Here N; is the total number of time steps. For simplicity, we consider only 
constant elements where collocation node r; is located at the center of the element and the nodal 
basis function is 


1, rs on element EF; that contains node r; 
yj(rs) = (42) 
0, otherwise 


The temporal basis function is taken to be the third-order shifted Lagrange basis polynomial that 
is commonly used for time domain boundary element methods:!* !® 


en(t) = 9 (| (43) 


where 
1+ ¢r+7°4+ 47° -l<7r<0 
1+ 57-77-3473 O0<Tr<l 
W(r)=4 1-gr—7? +579 Lae = (44) 
1-9r+r-i 2<7T<3 
0 other 
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For example, at a point rs on element E; and at an off-nodal time t = t, — nAt, 0 < n < 1, the 
value for $(1s,t) is found by 


(r,t) = girs) [wR (—n) + uF -NY(L = 9) + uF PW(2 — 9) + uF F(3 — 0) (45) 


With nodal spatial and temporal basis functions defined above, the expansion coefficient ue in (41) 
represents the value of ¢ at the collocation node rj; on element EF; at time level t,. By substituting 
expansion (41) into boundary integral equation (34) and evaluating the equation at collocation 
points r; of all elements and at time level t,, a March-On-in-Time scheme (MOT) is obtained that 
can be expressed in a matrix form as 


Bou” = q" — Byu""! — Bou”? —--- Byu”7 (46) 


where u* denotes a vector that contains all unknowns {ut ae Ne} at time level t,. The 


nonzero entries for matrices B;,, k = 0,1,2,..., J, in (46) are: 


: aie G5 ae: . ; 
{Bu}ij = 27007 Vp_1 (tn) + af Pa Cac) uy a) drs + bedi;0K0Di 
z Go : Bh oiis b =30G09Go jn an 
+be iA OnOn (von) — 6ij0n—K(tn) + a) drs + ae I. R an’ On n—k(tR)ars 


(47) 
for 7,7 = 1,2,...,Ne, where 6;; and 6,9 are Kronecker delta functions and a prime in the above 
denotes the derivative with respect to time and 


“= R(rs, 15 0G 
th =tn+B-(ri—Ts) ss ) D;=- L. Apion re riers (48) 


It is easy to see that the entry {B;}, ; Tepresents contributions to the value at node r; from the nodal 
value of element E; of time level t,_,;. The integrals in (47) are to be evaluated using high-order 
quadrature on each element. For the computational results reported in this paper, each element is 
mapped into a standard element of [—1, 1] x [—1, 1] and Legendre-Gauss quadrature rule of degree 
6 is used for integration in each dimension. Integration on the singular elements when 7 = 7 is 
detailed in Appendix B. 


The index J in (46) denotes the maximum time history of the solution required for (46) and is 
dependent on the length of the scattering surface and the mean flow as 


L = _ 
= noe AE +3, L= ae [—M : (r, _ rs) + R(rs,r,)| (49) 


Due to the limited temporal stencil width shown in (44) and (45), the B matrices are sparse. 
In particular, we note that matrix Bo in (46) is a very sparse matrix and represents interactions 
within the same element or between nearby nodes at the same time level t,. Bo is also found to be 
diagonally dominant. Solutions for w” in (46) can be found efficiently by an iterative method, such 
as the Jacobi iterative method, with rapid convergence.” !® 
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VII. Stability and eigenvalue analysis of the new integral equation 


As mentioned in previous sections, direct numerical solution of the time domain boundary integral 
equation (27) is prone to numerical instabilities. In Figure 2, we first show an example of scattering 
by a parabolic wing in a mean flow of Mach number 0.5, M = (0.5,0,0), to demonstrate the elim- 
ination of numerical instability by the Burton-Miller reformulation of TDBIE (27). The geometry 
of the scattering surface is a convex parabolic wing and is defined as follows: 


z=0.1L,(1—27/L2), -Lygsa<Ly, -ly<y<ly (50) 


where Lz = Ly = 0.5. The scattering surface is discretized by 2316 quadrilateral elements as shown 
in Figure 2. The source function is a broadband point source defined as the following: 


q(r,t) =e~% 5(r — r9) (51) 


where ro = (0,0, 1) and o = 1.42/(6At)?. The time history of the solution on a surface collocation 
point is plotted in Figure 2 for the cases without and with Burton-Miller reformulation. The top 
figure shows the result obtained by directly solving the TDBIE (27). It is seen that the solution 
behaves well initially but eventually becomes unstable. On the other hand, the solution obtained 
by the BM-TDBIE (34), shown in the bottom figure, remains stable. 


To further study the stability of the MOT scheme (46), we conduct a numerical eigenvalue study of 
the discretized system of equations.? We look for solutions of the form 


wu” = "en (52) 


to the corresponding homogeneous system of (46). By substituting (52) into (46), we obtain a 
polynomial eigenvalue problem 


[Bod\’ + By\’ 1 + Bod’? +--+ By-1\+ By] ep =0 (53) 


which can be cast into a generalized eigenvalue problem as follows: 


—-B, -Bo ::: -::. -Bj_1 —-By ej-1 Bo 0 0 tee 0 0 ey-1 
T 0 Rie Aho 0 0 ej—2 0 O. :: 0 0 es—2 
0 Fo 8s Gees 0 0 . -) 0 0 . 0 60 (54) 
0 0 anes ahs 0 0 e1 0 0 0 0 e1 
0 0 eats eats T 0 €0 0 0 O --. O I €0 


where e; = eo. Numerical scheme (46) is stable if |\| < 1 for all eigenvalues of (54). 


Eigenvalue analyses of scattering by two geometric shapes are presented in Table 1. One of the 
geometries is the parabolic wing as described previously in (50). The other is a sphere of radius 
0.5. The surface of the sphere is first discretized by 512 unstructured triangular elements each of 
which is then subdivided into three quadrilateral surface elements resulting in a total of 1536 surface 
elements. The mean flow Mach number varies from 0 to 0.9. A total of eight cases are considered 
in Table 1. 


Eigenvalues of the generalized eigenvalue problem (54) can be found via a sparse eigenvalue solver 
available in MATLAB and Python, or by a matrix power iteration method detailed in Appendix C. 
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eae Without Burton-Miller reformulation 
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Figure 2. Time history of numerical solution on a surface collocation point, showing the elimination 
of instability by Burton-Miller reformulation of TDBIE. M = (0.5,0,0), At = 0.02. Top: solution of 
(27) without Burton-Miller reformulation; bottom: solution by BM-TDBIE (34). 
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The values of the largest eigenvalue for the eight cases are listed in Table 1. For the Burton-Miller 
formulation BM-TDBIE (34), all eigenvalues are no greater than unity and stability is observed. In 
contrast, direct solution of the equation (27) results in eigenvalues greater than unity in all but two 
of the eight cases studied, indicating that equation (27) without Burton-Miller reformulation can 
lead to unstable solutions. 


Table 1. Maximum eigenvalue, |A|,,,,.. computed by (54) for scattering by a parabolic wing and a 
sphere, for cases with and without Burton-Miller (B-M) reformulation. N. is the total number of 
elements and M is the mean flow Mach number. 


Parabolic Wing Sphere 
Minas Minas 

Ne M with B-M without B-M Ne M with B-M without B-M 

Eq. (34) Eq. (27) Eq. (34) Eq. (27) 
2316 | 0.0 1.000000 1.095949 1536 | 0.0 1.000000 1.007840 
2316 | 0.3 1.000000 1.160628 1536 | 0.3 1.000000 1.000000 
2316 | 0.6 1.000000 1.129116 1536 | 0.6 1.000000 0.999968 
2316 | 0.9 1.000000 1.582909 1536 | 0.9 1.000000 1.003901 


VIII. Conclusions 


In this paper, we have considered the boundary condition to be used in time domain boundary inte- 
gral equation analysis of acoustic scattering by solid bodies under a constant mean flow assumption. 
After an examination of the energy equation associated with the convective wave equation, it is pro- 
posed that the boundary condition be defined by the requirement that the energy flux be zero at 
solid boundaries. A new TDBIE is derived based on this new solid wall boundary condition. The 
new formulation differs from those found in the literature on the part of the boundary where the 
constant mean flow itself does not satisfy the solid surface boundary condition. In addition to 
conserving the acoustic energy, another significant advantage of the new equation is that it is con- 
siderably simpler than previous formulations. In particular, tangential derivatives of the solution on 
the solid surfaces are no longer needed in the new formulation, which greatly simplifies numerical 
implementation. To stabilize the TDBIE, Burton-Miller reformulation is also derived. Numerical 
solutions and eigenvalue analysis are presented that demonstrate stability of the new formulation. 


Naturally, from a physical point of view, the null acoustic energy flux condition at rigid surfaces 
should be equivalent to, or a direct consequence of, the condition that the normal acoustic velocity 
becomes zero on rigid surfaces. The fact that the two now differ in the formulation of the boundary 
integral equation for scattering with flow is due to the inconsistency on the part of the underlying 
mean flow itself when the constant flow simplification is made. Thus, boundary integral equation 
approaches with a constant mean flow would be applicable only to problems where such a sim- 
plification is acceptable or justified, such as in scattering with flow over slender bodies. From a 
computational point of view, however, the current formulation based on the energy flux condition is 
significantly simpler than those based on the normal velocity condition. As such, as a result of the 
present analysis, the enforcement of normal acoustic velocity being zero on boundary points where 
the mean flow itself does not satisfy such a condition, and the computational complications it brings 
in with the separation of normal and tangential derivatives of the solution, become unnecessary. 
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Appendix 


A. Limit of weakly-singular integral 


By equations (17) and (32), it is easy show that the modified normal derivative Go (rs,r,) and 


s 
Go (r;,7,) have a singularity of order O(1/|r; — r4,|) and O(1/|r, — r,)|°), respectively, which 
makes their surface integrals weakly-singular and hyper-singular, respectively. In this appendix, we 


state some useful results. 


For surface integrals involving 5 Go , we have 


r’' €V, exterior of 9 
= 7CES (55) 
r’€V_, interior of S 


OGo 


¢ OF —(r,,r')dr, = 


ewe © 


The first and third equations in (55) can be obtained by the fact that any constant can be a solution 
to the homogeneous convective wave equation with homogeneous normal derivative on the boundary 
for the interior domain V~ enclosed by S. By substituting p = 1 to equation (15) and noting the 
choice of the normal direction and the placement of r’, the first and third equation in (55) follow 
immediately. 

The second integral in (55) becomes weakly singular when r’ approaches a point on surface $. This 
particular limit has been studied previous in the literature for a mean flow that is aligned with the 
x-coordinate.?!;?6 Here, we show the calculation for a general mean flow. Assuming r’, is a smooth 
point on S, consider modifying surface S by a spherical surface of radius € and centered at r‘, as 
shown in Figure 3. The surface is assumed to be smooth at r/,. If we denote the small hemispherical 
surface as S., we have 


lim CON gcd: = lim ED ogee. lim 260 (n,.r!)dr (56) 
n 


rior, Jg On rr, JS—S. On rr, JS. 


Note that, for the surface integral on S,, using (10), we have 


0Go __amiles— sh) + nolys— 9h) +a(2— 24) __ oa € 


On RS a? 
By the symmetry of R with respect to hemispheres S$, and $%, the complementary hemisphere of S_, 
and by using a local spherical coordinate system centered at r/, and its local z direction coincides 
with mean flow M , namely x, — 2, = esinvcos9, y, — y, = esinv sin, z, — 2, = ecosv, we have 


fF Go , 2 € d a? / € Ta e sin 
im —dr,=-—a =.dr, = —-— = =-—_ 

3 gd 3/2 
rir, Jg. On elt 2 Jse+S! R €2 cos? v + €2a2 sin? v) / 
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duvdé 


1 
1 
-1 (a? + (1 — a?)x2)°? 


The last integral above can be found by direct integration. The second equation in (55) follows as 
€ — 0 and by noting that, for r’ € V, the limit on the left hand side of (56) is zero. 


Figure 3. A schematic diagram for a hemisphere that caps a surface point r’,. Note that the normal 
vector is in the direction outward from the region of solution and into the body. 


B. Evaluation of hyper-singular integral 
We consider the numerical evaluation of the regularized integral involving the double normal deriva- 


tive of Go in (34). Note that as r, > ri, 


/ / / R 0 / / / / / 0 / / / 
(rast) Ora tay SE (asthe) = VOC at) (me) +B) so (Mla) +O( re w4?) (87) 


Let a surface element E; be mapped to a local coordinate (€,7) € [—1, 1] x [-1, 1], which is then in 
turn converted into a local polar coordinate (r,@) centered at the collocation point r/,. Denote the 
integrand for the integral in (r,@) as 


0G ‘ RO : 
F(r8) = (SEE) (orate) — oleat) + 2s (rate) bee % to (58) 
By (57), F(r,0) is of order O(1/r) as r > 0. Let 
lim r?F(r,0) = G(0) (59) 


It is easy to show that fa G(0)d@ = 0. Then we have the following for the integral on surface 
element Ej: 


Qn Qxr pr(O) 2 _ 
lim i a (r,0)rdrd@ = lim ie, ie eee) | OO ids 
r r 


e>0 e>0 


Qa (9) » 2 Qn 
=| i HG) 64 rlim | G(0)[Inr(0) — Ineldé 
E> 0 


Q2r pr(O) 12 _ Qr 
: / / EMG 2) CO arao + [  G(0)nr(6)d0 
0 


The final integrals can be evaluated using high-order numerical quadrature. 
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C. Eigenvalue by matrix power iteration method 


We describe a matrix power iteration method for finding the largest eigenvalue of (54). Let 


—B)'B, —Bj'B2 --. «+ —Bo'Bs-1 —Bo'Bs 
T 0 ahd abe 0 0 
Gee | 0 T ik Asay) 0 0 | 
0 0 eevee 0 0 
0 0 iat iad T 0 
Then, the power iteration method proceeds as follows:? 
given an arbitrary unit vector e), and for k = 1,2,..., compute 
vo) = Ae*-) (60) 
(k) 
ian © 
ee) = 61 
[eM vs 
and eigenvalue 
A) = [e®] * Ae) = [e] P ay(h+) (62) 


The iteration is stopped when |) — NERD rs |A®| < ¢€, where € is the tolerance and set to be 
10-!?. When the iteration is convergent, it converges to the largest eigenvalue of A. 


In particular, if we denote 


€ 7-2 J—2 

efh) — : , Vl) = , (63) 
eff) 7 (t) 
elt) ‘ (t) 


then, equation (60) can also be computed through the following relations that save memory and 
storage: 


Ca —B;* [Bie + Boe,” feeet Bye?) + Bye?) 0), = Eee vee uh) — ek) 
(64) 


We note that the iterative step shown in (64) is the same as the MOT iteration (46) without the 
source term. Therefore, it can also carried out using the same computational scheme for (46). 


References 


TA. J. Burton and G. F. Miller, The application of integral equation methods to the numerical solution of some 
exterior boundary-value problems, Proc.R.Soc.London, Ser A, Vol. 323, 201-210, 1971. 


18 of 19 


American Institute of Aeronautics and Astronautics 


2D. J. Chappell, P. J. Harris, D. Henwood and R. Chakrabarti, “A stable boundary element method for modeling 
transient acoustic radiation”, Journal of Acoustical Society of America, Vol. 120, No. 1, 74-80, 2006. 


3S. J. Dodson, 8. P. Walker and M. J. Bluck, Impicitness and stability of time domain integral equation scattering 
analysis, The Applied Computational Eletromagnetics Society J., Vol. 13, 291-301, 1998. 

“Dowling, A. P. and Ffowcs Williams, J. E., Sound and Sources of Sound, Horwood Publishing, Westergate, 1983. 

5M. H. Dunn and A. F. Tenetti, “Application of fast multipole methods to the NASA Fast Scattering code”, AIAA 
paper 2008-2875, 2008 


SA. A. Ergin, B. Shanker and E. Michielssen, “Analysis of transient wave scattering from rigid bodies using a 
Burton-Miller approach”, Journal of Acoustical Society of America, Vol. 106, No. 5, 2396-2404, 1999. 

"P. Farassat and M. K. Myers, Extension of Kirchhoff formula to radiation from moving surfaces, Journal of Sound 
and Vibration, Vol. 123, 451-460, 1988. 


8J. E. Ffowcs Williams and D. L. Hawkings, Sound generation by turbulence and surfaces in arbitrary motion, Phil. 
Trans. Royal Soc., Vol. 264A, 321-342, 1969. 

°G. H. Golub and C. F. Van Loan, Matrix Computation, 4th Edition, Johns Hopkins Studies in the Mathematical 
Sciences, 2013. 

10Y Guo, Computation of sound propagation by boundary element method, NASA Contract Report, NAS1-00086- 
A003, 2005. 

UF. Q. Hu, “A spectral boundary integral equation method for the 2D Helmholtz equation”, Journal of Computational 
Physics, Vol. 120, 340-347, 1995. 

2A Q. Hu, “A fast numerical solution of scattering by a cylinder: spectral method for the boundary integral 
equations”, Journal of Acoustic Society of America, Vol. 96, 3639-3703, 1994. 

13% Q. Hu, An efficient solution of time domain boundary integral equations for acoustic scattering and its acceler- 
ation by Graphics Processing Units , AIAA paper 2013-2018, 2013 

MF Q. Hu, Further Development of a Time Domain Boundary Integral Equation Method for Aeroacoustic Scattering 
Computation , AIAA paper 2014-3194, 2014. 

15% Q. Hu, Y. P. Guo and A. D. Jones, On the computation and application of exact Green’s function in acoustic 
analogy , AIAA paper 2005-2986, 2005. 

16F Q. Hu, M. E. Pizzo, and D. M. Nark, On the assessment of acoustic scattering and shielding by time domain 
boundary integral equation solutions , AIAA paper 2016-2779, 2016. 


'TH-W Jiang and J-G Ih, Stabilization of time domain acoustic boundary element method for the exterior problem 
avoiding the nonuniqueness, Journal of Acoust. Soc. Am, Vol. 133, 1237-1244, 2013. 

18A.D. Jones and F. Q. Hu, “A three-dimensional time-domain boundary element method for the computation of 
exact Green’s functions in acoustic analogy”, AIAA paper 2007-3479, 2007. 

'Sy_W. Lee and D. J. Lee, “Derivation and implementation of the boundary integral formula for the convective 
acoustic wave equation in time domain”, Journal of the Acoustical Society of America, Vol. 136, 2959-2967, 2014. 
201. P. Lockard, An efficient, two-dimensional implementation of the Ffowcs Williams and Hawkings equation, 
Journal of Sound and Vibration, Vol. 229, 897-911, 2000. 

217, Long, The compressible aerodynamics of rotating blades based on an acoustic formulation, NASA TP 2197, 
1983. 


22\W. Mohring, Energy conservation, time reversal invariance and reciprocity in ducts with flow, Journal of Fluid 
Mech., Vol. 431, 223-237, 2001. 


23\V. R. Morgans, The Kirchhoff formula extended to a moving surface, Philol. Mag., Vol. 9, 141-161, 1930. 
240. L. Morfey, Acoustic energy in non-uniform flows, Journal of Sound and Vibration, Vol. 14, 159-179, 1971. 
25P. M. Morse and K. U. Ingard, Theoretical Acoustics, Princeton, 1986. 


26. K. Myers and J. S. Hausmann, On the application of the Kirchhoff formula for moving surfaces, Journal of 
Sound and Vibration, Vol. 139, 174-178, 1990. 

27M. K. Myers, Transport of energy by disturbances in arbitrary flows, Journal of Fluid Mechanics, Vol. 226, 383-400, 
1991. 


°8M. K. Myers and J. S. Hausmann, Computation of acoustic scattering from a moving rigid surface, Vol. 91, 2594- 
2605, 1992. 

29H. A. Schenck: Improved Integral Formulation for Acoustic Radiation Problems, Journal of Acoustical Society of 
America , Vol .44, 41-58, 1968. 

3°P, Zhang and T. W. Wu, A hypersingular integral formulation for acoustic radiation in moving flows, Journal of 
Sound and Vibration, Vol. 206, 309-326, 1997. 


19 of 19 


American Institute of Aeronautics and Astronautics 


